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ABSTRACT 

We add black holes to nonrotating, spherical galaxy models, with the assumption 
that the black-hole growth is slow compared with the dynamical time but fast compared 
with the relaxation time. The outcome differs depending on whether the core of the 
initial galaxy does or does not resemble that of an isothermal sphere. For the isothermal 
case the previously-known results are confirmed and sharpened: the black hole induces 
cusps in the density (p ~ r~"^'^) and velocity dispersion (v"^ ~ r~^), and a tangential 
anisotropy in the velocity distribution away from the center. For the non-isothermal case 
the induced density cusp is steeper, and the induced anisotropy is larger and penetrates 
right to the center. The cusp around the black hole is insensitive to anisotropy in the 
initial velocity distribution, and also to the origin of the black hole, unless its mass 
comes exclusively from the stars of lowest angular momentum, in which case the cusp 
is suppressed. We discuss the implications for the interpretation of evidence for massive 
black holes in galactic nuclei. 



1. Introduction 

The growth of central mass concentrations appears 
to be a natural result of the processes that shape large 
stellar systems. Quasars and active galactic nuclei are 
believed to derive their power from the most spec- 
tacular outcome of this: matter at the center of a 
galaxy that has collapsed into a massive black hole 
(BH). If this view is correct then many galaxies should 
today contain "dead quasars," massive BHs starved 
of fuel (Rees 1990). The search for dynamical evi- 
dence of BHs from the distribution and kinematics 
of stars in the centers of galaxies has thus received 
considerable attention, as constraints on BH masses 
would help us understand the structure of quasars 
and, more generally, the formation and evolution of 
dense galactic nuclei. There is evidence from ground- 
based observations that BHs have been detected in 
this way in a half-dozen or so nearby galaxies (see 
the reviews by Dressier 1989, Gerhard 1992, and Kor- 
mendy 1992), although a watertight case has yet to be 
made. Some of the remaining questions might be an- 
swered by high-resolution observations with the refur- 
bished Hubble Space Telescope (HST); others could 
perhaps be answered by better theoretical modelling 
of existing data. In this paper we avoid detailed mod- 
elling of particular galaxies, and focus instead on a 
complementary approach — the construction of gen- 
eral theoretical models of galaxies with central BHs 
— in the hope that this might help both in identify- 
ing and measuring BHs in galaxies, and in explaining 
the origin of these systems. 

By the construction of models of galaxies with cen- 
tral BHs we mean the exploration of equilibrium solu- 
tions that result from plausible initial conditions and 
definite (though perhaps speculative) BH formation 
scenarios. We thus exclude models in which a BH is 
placed at the center of a stellar system without regard 
to the origin of that configuration, with simplifying 
assumptions made without justification. Examples 
of these include the "loaded polytropes" of Huntley 
and Saslaw (1975), which assume an isotropic veloc- 
ity distribution and a polytropic relation between the 
pressure and density, and the "ry-models with BHs" of 
Tremaine et al. (1994), which assume an isotropic ve- 
locity distribution and a law for the variation of den- 
sity with radius. These models can give helpful math- 
ematical insights into the range of allowed solutions, 
but their arbitrary nature makes it unlikely that they 
will match real galaxies well. We also exclude models 



constructed by techniques such as linear programming 
to match observations of galaxies believed to contain 
BHs (see, e.g., Richstone and Tremaine 1985). These 
models are valuable for deciding whether particular 
galaxies contain BHs (and for constraining the BH 
masses if they do), but they are not based on initial 
conditions and a BH formation scenario, and hence 
can neither help us assess theories for these nor make 
general predictions for the solutions we should expect 
to find. 

The expected distribution of stars around a mas- 
sive BH was studied in detail in the late 1970s and 
early 1980s, but in the context of globular clusters, 
not galactic nuclei (see Shapiro 1985 for a review). 
The key assumption underlying this work is that the 
two-body relaxation time is short compared with the 
age of the system. The cluster is then driven by re- 
laxation processes to a steady-state solution in which 
the consumption of stars by the BH and the diffusion 
of new stars into the loss cone result in a density cusp 
p(r) ~ r-^/4 (Bahcall and Wolf 1976). 

In most galaxy cores, however, the relaxation time 
is long compared with the age of the system, and there 
is no reason to expect a unique solution. The distribu- 
tion of stars around a central BH in a galaxy is likely 
to depend upon many factors, including the order in 
which the galaxy and the BH form, the structure of 
the galaxy — assuming that it formed first — before 
the BH forms (spherical, axisymmetric, or triaxial? 
rotating or nonrotating? cuspy or fiat? isotropic or 
anisotropic?), and the origin of the BH mass (stars? 
gas? an external source?). If we had predictions 
for the cusps that result from all the possibilities we 
would have a good understanding of what we might 
find, and could rule some possibilities out by compar- 
ing their predictions with observations. 

Only one plausible formation scenario has been ex- 
plored in detail: the growth of a central BH from 
the accumulation of gas on a time scale long enough 
that the stellar action variables are adiabatically con- 
served. The consequences of this can be derived most 
easily for spherical galaxies. Peebles (1972) consid- 
ered the adiabatic growth of a central BH in an 
isothermal sphere, and showed it would lead to a den- 
sity cusp p(r) ~ r~^". Young (1980) constructed 
numerical models that confirmed Peebles's result and 
showed that the BH induces a tangential anisotropy 
in the velocity distribution. Goodman and Binney 
(1984) obtained an approximate solution to the same 
problem, and showed that the velocity distribution 



remains isotropic at the center despite the tangential 
bias induced nearby. Lee and Goodman (1989) gen- 
eralized Young's calculation in an approximate way 
to study the growth of BHs in axisymmetric, rotating 
galaxies. 

The adiabatic growth of a BH in a triaxial galaxy 
is more difficult to analyse than the corresponding 
problem for a spherical galaxy, but is potentially more 
interesting because the maintenance of triaxiality de- 
pends on the existence and selective population of 
box orbits, orbits that — given enough time — pass 
arbitrarily close to the central BH. Most studies of 
this problem have considered only the effect of a BH 
on single-particle orbits, and have not attempted to 
follow the self-consistent evolution of a whole galaxy. 
Gerhard and Binney (1985) argued that the scatter- 
ing and redistribution of box orbits by a BH would 
force the inner regions of a triaxial galaxy to become 
rounder, out to a distance of about 1 kpc for a 10® Mq 
BH in a giant elliptical galaxy, perhaps resulting in 
a global change in the galaxy shape. Similar argu- 
ments were made by Hasan and Norman (1990) for 
the growth of BHs in barred galaxies. Norman, May, 
and van Albada (1985) pioneered the use of N-body 
simulations to follow the self-consistent evolution of 
triaxial galaxies with growing BHs, and found results 
consistent with the predictions of Gerhard and Bin- 
ney. The number of particles in their simulations was 
small (5000-20000), however, which limited the real- 
ism with which they could model the systems, and 
which caused spurious relaxation that influenced the 
results in an uncertain way. 

Despite these advances, many questions remain 
about the effect a growing central BH has on the 
structure of a galaxy. This is true even for spherical 
galaxies. We don't know what extremes are possible, 
e.g., how steep or gradual the density and velocity 
cusps can be, and whether it is possible to hide a BH 
at the center of a galaxy without an observable cusp. 
Nobody has made accurate predictions for the tangen- 
tial anisotropy expected near a BH, which we need to 
interpret observations and to make precise mass esti- 
mates. We don't know exactly what signature a BH 
will leave in the distribution of velocities along the 
line-of-sight, a pressing question now that an effort 
is underway to quantify the deviation of these dis- 
tributions from Gaussians and to use that informa- 
tion to constrain dynamical models (van der Marel 
et al. 1994). Another question is to what extent a 
BH can suppress instabilities such as the radial-orbit 



instability that plague some models for galactic nu- 
clei without BHs (see, e.g., Merritt 1987; Palmer and 
Papaloizou 1988). 

For non-spherical galaxies the questions are more 
numerous and profound. The main question is whether 
a central BH is compatible with triaxiality, i.e., whether 
a the growth of a BH in a triaxial galaxy forces the 
galaxy to become rounder and, if it does, whether this 
happens gradually from the inside out, or abruptly in 
a manner that affects the whole galaxy. The answer 
might allow us to constrain BH masses and formation 
histories from observed isophote shapes, and might 
help explain why unresolved elliptical galaxy cores 
tend to appear disky while resolved ones tend to ap- 
pear boxy (Nieto, Bender, and Surma 1991). Figure 
rotation and resonances can mitigate the influence of 
the BH and complicate the analysis by preventing or- 
bits from approaching close to the center (Pfenniger 
and de Zeeuw 1989). The growth of a central BH in a 
triaxial galaxy will cause many orbits that were regu- 
lar to become stochastic, although it is not clear what 
stochasticity implies for the structure of the galaxy 
(Udry and Pfenniger 1988). Many of these questions 
are difficult and can be addressed only with the help 
of large N-body experiments. 

In this paper we start with the simplest problem, 
and study the adiabatic growth of a central BH in a 
nonrotating, spherical galaxy using the numerical ap- 
proach of Young (1980). There are several reasons for 
re-examining this problem. The first is that previous 
work considered the growth of a BH in only one galaxy 
model, the isothermal sphere, and it is not clear which 
of the results are peculiar to this model and which 
are general. We now know that few elliptical-galaxy 
cores resemble that of the isothermal sphere; many 
have surface brightnesses that continue to rise at the 
smallest observable radii. We examine a family of 
simple galaxy models with this property, and show 
that the adiabatic growth of a central BH gives re- 
sults for them that differ qualitatively from those for 
the isothermal sphere. In particular, the density cusp 
induced by the BH is steeper, and the anisotropy in 
the velocity distribution is larger and penetrates right 
to the center. 

Another reason for re-examining this problem is 
that previous work did not consider the origin of the 
BH, and assumed either that its mass came from a 
source external to the stars (such as gas that seeps in 
from the outer parts of the galaxy), or, if the BH did 
grow at the expense of the stars, that the mass loss 



could be ignored. We perform calculations both with 
and without taking stellar mass loss into account, and 
show that the mass loss is indeed unimportant unless 
it is highly concentrated towards the galaxy center. 

Finally, we are re-examining this problem to ex- 
tract quantitative results that were ignored or con- 
sidered only qualitatively in previous work, such as 
the the anisotropy in the velocity dispersion, and the 
fourth moment of the line-of-sight velocity distribu- 
tion. We are doing this partly because of their im- 
portance for the interpretation of observations, but 
also because we want to use them to calibrate an N- 
body program we are developing to study the growth 
of BHs in triaxial galaxies, a problem that cannot be 
handled by the simple techniques used here. 

We assume that the mass concentration at the 
galaxy center is a massive BH, although some of our 
conclusions will apply to galaxies with other mass 
concentrations, such as star clusters, provided that 
they are sufficiently dense. 

2. Computational methods 
2.1. Strategy 

The computational strategy and equations on which 
our calculations are based are described in detail by 
Young (1980), and are summarized only briefly here. 

The starting point of the calculation is an equilib- 
rium spherical galaxy with no BH (or a small BH). 
Young then adds a BH to the center making two as- 
sumptions: flrst, that the BH mass comes from a 
source external to the stellar system, and does not 
deplete the stellar distribution function; second, that 
the BH growth is slow enough that the stellar action 
variables (radial action J^ and angular momentum L) 
are adiabatically conserved, but fast enough that two- 
body relaxation can be ignored. The flrst assumption 
is not essential and can be relaxed; we do this later 
in the paper and flnd that the results are not highly 
sensitive to the origin of the BH mass. The second as- 
sumption is necessary in this type of calculation, but 
is justifled provided that the relaxation time is long 
compared with the age of the system, and that the 
BH growth is slow compared with the orbital period 
near the galaxy center. 

Young's approach does not follow the time-evolution 
of the system, but solves directly for the flnal self- 
consistent distribution of stars in which the distribu- 
tion function has the same dependence on the action 



variables as it did in the original model. This is done 
by adding a central BH and then passing repeatedly 
through the following loop to converge onto the so- 
lution: compute the potential from the mass distri- 
bution (including the BH); compute the action vari- 
ables in that potential; adjust the distribution func- 
tion f(E, L) to remain a flxed function of the actions; 
compute the density generated by the distribution 
function in the current potential; check how much the 
density has changed from the previous iteration; de- 
cide whether to accept the solution or to pass through 
the loop again. A dozen or so repetitions are usually 
enough to reduce the change in density to less than 
one part in 10'* at all radii, which is our convergence 
criterion. 

2.2. Computer program and output 

Our program that implements this strategy is sim- 
ple and easy to run, yet flexible enough to handle a 
variety of galaxy models. The density, distribution 
function, and other properties of the galaxy are de- 
scribed by their values on a discrete set of grid points: 
Pi = p(ri), fij = f(Ei,Lj), etc. The radial grid 
points are spaced logarithmically between minimum 
and maximum values chosen by the user (typically 
Vynin = 10"'* and fmax = 10^ in our units). The 
grid points for energy are chosen to match the po- 
tential at the radial grid points, and thus vary during 
the iterative procedure to converge onto the new po- 
tential. The grid points for angular-momentum are 



spaced linearly between x„ 



and Xr, 



1, 



where x = L/Lc is the ratio of the angular momen- 
tum to the circular angular momentum at the given 
energy. For the calculations in this paper we used 200 
grid points for radius and energy and 20 for angular- 
momentum, although that is more than is necessary: 
the gross properties of the stellar cusps can be repro- 
duced easily with half this number. A typical calcula- 
tion takes about one minute to complete on our IBM 
580 RISC workstation. 

At the end of the calculation we have the complete 
distribution function for the stellar system, which we 
condense into a small number of moments as a func- 
tion of radius (r) or projected radius (R). We flrst 
compute the density 
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where Lm is the maximum angular momentum at- 
tainable by an orbit of energy E at radius r, 
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and Vr and Vt are the radial and tangential velocities, 
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We then project the intrinsic density and velocity mo- 
ments onto the plane of the sky to get the surface 
density 
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the dispersion of the line-of-sight velocity distribution 
(LOSVD) 
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and the fourth moment of the LOSVD (Merrifield and 
Kent 1990) 
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which we usually present as the dimensionless kurto- 
sis K = {Vp)/{Vp)^- We also compute the anisotropy 
parameter 

..2\ //r)„.2\ 



P=l-{vi)/{2vi), 



(8) 



which is for an isotropic distribution and can vary 
between — cxd (for purely circular orbits) and -|-1 (for 
purely radial orbits). 

Some of our output quantities can be compared 
directly with observations and some cannot. The 
anisotropy parameter /3 cannot, but it is nevertheless 
important because it must be known (or assumed) be- 
fore observed values for S and (v?,) can be converted 
into a precise mass estimate. The velocity moments 
(Vp) and (Vp) could be compared directly with obser- 
vations if we could observe without noise and with 



infinite resolution, but in practice we cannot and the 
moments are affected by noise in the wings of the dis- 
tribution (this is especially true for (Vp), which can- 
not be measured reliably) and by the nonzero see- 
ing radius over which the observations are averaged. 
Gerhard (1993) and van der Marel and Franx (1993) 
argue that it is better to quantify observed LOSVDs 
by a set of Gauss-Hermite moments, which describe 
the LOSVD by a Gaussian fit and the deviations from 
this fit, and which are not as sensitive to the wings of 
the distribution as are moments such as (Vp) and (Vp). 
It might have been interesting to present our output 
in this form, but that would have made the program 
more complicated. The classical moments (Vp) and 
(Vp) are easy for us to compute and are sufficient to 
give an understanding of the intrinsic dynamics. For 
distributions close to a Gaussian, the Gauss-Hermite 
moment /}4 of van der Marel and Franx (1993) is re- 
lated to the kurtosis by 

K ~ 3 -F 8^6/14 (9) 

(this linear relation is unreliable if |/}4| ^ 0.03). 

2.3. Tests of program 

To test the program we first ran it with simple 
galaxy models (described in Section 3.1.) with no BH 
to check that the output quantities matched those of 
the input models to sufficient accuracy (at least one 
part in lO'*, except at the outermost radial grid points 
where the accuracy is worse) and that the accuracy 
improved in the expected manner when the number of 
grid points was doubled or quadrupled. To check the 
calculations of the anisotropy parameter /3 and kur- 
tosis K we used the anisotropic Plummer models of 
Dejonghe (1987) and Cuddeford (1991), and several 
other anisotropic models that we derived by Cudde- 
ford 's technique. 

The main test of the program was to reproduce 
Young's (1980) results for the adiabatic growth of a 
BH at the center of an isothermal sphere. A visual 
comparison of our results with his showed satisfactory 
agreement. The only quantity for which we could 
detect any disagreement is the ratio /max//min = 
f(E, Lc(E))/f(E, 0), which Young plotted in his Fig- 
ure 1 for just one BH mass (the largest he considered) 
to quantify the anisotropy in the distribution func- 
tion. We show our version of this plot in Figure 1(a), 
in the same units used by Young. Our /max//min ratio 
approaches unity towards the left of the plot slightly 
slower than Young's does, but the difference is small 



and we do not view it as significant (our result does 
not cliange if we double or quadruple the grid resolu- 
tion). 

It is unfortunate that Young presented the anisotropy 
in this way, as it is not obvious how large an anisotropy 
it implies for the velocity dispersion (although Young 
concluded correctly that even for his largest BH mass 
the anisotropy was small). We show in Figure 1(b) 
the anisotropy parameter (3 for the same calculation 
as in Figure 1(a). Note first that the anisotropy is 
small, as Young concluded, and second that it goes 
to zero at the center, as predicted by Goodman and 
Binney (1984). The Goodman-Binney solution^ gives 
an anisotropy that varies with radius as /3 ~ r^'"^, 
which we show in Figure 1(b). This solution is ex- 
pected to be valid only near the galaxy center, as it 
is derived by approximating the true potential by a 
harmonic potential before the BH is added and by a 
Kepler potential after. Binney and Petit (1989) es- 
timate the boundary to the solution's validity to be 
the radius at which the enclosed stellar mass in the 
initial model equals the mass of the BH, but for the 
calculation shown in Figure 1 we find that at this ra- 
dius (r ~ 0.6) the anisotropy implied by the solution 
is about five times too large. 

The dotted line in Figure 1(b) shows our program 
output when we turn off the self-consistent potential 
calculation and impose harmonic and Kepler poten- 
tials at all radii (chosen to match the central poten- 
tials before and after the BH is added). The output 
agrees with the Goodman-Binney solution, confirm- 
ing that our program is correctly keeping the distribu- 
tion function a fixed function of the action variables. 

3. Results 

3.1. Initial models and output figures 

As the starting point for our calculations we pick 
well-known, simple galaxy models with mass distribu- 
tions that span the range of behaviours expected for 
spherical galaxies. The models are not intended to 
be accurate representations of real galaxies, although 
two of them (the 7 = 1 and 7 = 2 models) do give 
reasonable fits to an i?^''* law. 

The first model we pick to have a core like that 



There is a typographical error in the Goodman-Binney paper 
suggesting that /? varies as r and not r ' . The Xm in their 
equation (13b) should be x'^. The same error is repeated by 
Binney and Petit (1989). 



of the isothermal sphere. We do not use the isother- 
mal sphere, because it was studied in detail by Young 
(1980). Instead we use Henon's (1960) isochrone 
model, defined by the potential 



(t>i{r) 



GM 
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The density corresponding to this potential falls off at 
large radii as r~'*. At small radii the density is nearly 
constant, and can be expanded in even powers of r: 



Pir) ~ KO) + \p"{^y + 



(11) 



Models that share this property — such as the Plum- 
mer model. King models, and the isothermal sphere 
— are often called models with isothermal cores, a 
name we dislike because the word isothermal should 
refer to the velocity distribution, not the density, and 
because the singular isothermal sphere has a steep 
density cusp {p ~ r~^) and yet certainly deserves to 
be called isothermal. It is nevertheless useful to have 
a name for galaxy models with the property (11), be- 
cause, as we shall soon see, they all respond in a sim- 
ilar manner to the adiabatic growth of a central BH. 
We shall call them models with analytic cores^, be- 
cause a density with spherical symmetry must be ex- 
pandable about the center as in (11) if it is an analytic 
function of the three spatial coordinates. 

For models of galaxies with non-analytic cores we 
pick three from the one-parameter family studied by 
Dehnen (1993) and Tremaine et al. (1994). We call 
them "7 models" because they are defined by the den- 
sity which Dehnen writes as 
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(Tremaine et al. describe the same density by the pa- 
rameter ry = 3 — 7 and call the models "ry models".) 
At large radii the density falls off as p-y(r) ~ r~'*, just 
as for the isochrone model, but at small radii the den- 
sity has a cusp p-y(r) ~ r~'* . The mass distribution 
Mj(r) is nonsingular as long as 7 < 3: 



M^(r) = M 



r + a 



3-7 



(13) 



Two of these models are well known from previous 
work: the 7=1 model is the Hernquist (1990) model; 



Suggested to us by S. Tremaine. 



7 = 2 is the Jaffe (1983) model. We use these as 
representative models for galaxies with mild (7 = 1) 
and steep (7 = 2) density cusps, and add as a third 
the 7 = model which, though it has a finite central 
density, does not qualify as a model with an analytic 
core because its density varies linearly with radius 
near the center. We have experimented with the 7 = 
3/2 model, but do not show the results here because 
they are intermediate between those for 7 = 1 and 
7 = 2 and do not reveal any surprises. 

The results from our calculations are shown in Fig- 
ures 2-5. We present all the results (except those 
for the isothermal sphere in Fig. 1) in the standard- 
ized units of Heggie and Mathieu (1986), in which the 
gravitational constant G and the initial galaxy mass 
M are both chosen to be unity, and the initial energy 
is chosen to be £" = —1/4. The scale lengths for the 
7 and isochrone models are thus a = 1/(5 — 27) and 
h = (37r — 8)/6. Each figure has four panels, showing 
(a) the surface density, (b) the projected velocity dis- 
persion, (c) the anisotropy parameter /3, and (d) the 
kurtosis of the LOSVD. Each panel of each figure has 
six lines: the dotted line shows the initial model be- 
fore the BH is added; the five solid lines show the final 
models after the adiabatic growth of BHs of masses 
0.001, 0.003, 0.01, 0.03, and 0.1. 

3.2. Surface-density cusps 

The surface-density cusps shown in Figures 2-5 
vary from one model to another, and, for all mod- 
els but the isochrone, are steeper than the S ~ R~^'''^ 
cusp found by Peebles (1972) and Young (1980) for 
the isothermal sphere. Varying the mass of the cen- 
tral BH merely shifts the radius where the limiting 
power-law cusp appears. 

The fact that the 7 models develop cusps steeper 
than S ~ i?~^" is perhaps not surprising for the 
models with 7 > 0, which have density cusps be- 
fore the BH is added, but it is for the 7 = model, 
which starts with a finite central density as does the 
isochrone and yet develops a cusp that is twice as 
steep. There is a simple explanation for this differ- 
ence. In deriving his cusp formula, Peebles (1972) 
assumed that the distribution function could be ap- 
proximated by a constant in the core of the initial 
model. This assumption is not valid for the 7 mod- 
els, because their distribution functions diverge as E 
approaches (t)(0). It is easy to generalize the deriva- 
tion to take this into account (see Appendix A). We 
need just three assumptions: that the initial model 



TABLE 1 

Adiabatic Density Cusps 



Model 



7 



A 



C 



isochrone 3/2 9/4 

7 = 00 1 2 9/4 

7=1 1 5/2 7/3 7/3 

7 = 3/2 3/2 9/2 12/5 12/5 

7 = 2 2 — 5/2 5/2 



has an isotropic core; that the potential varies with 
radius near r = as a power law (j> ~ r^"'''; and that 
the distribution function diverges near E = (t)(0) as a 
power-law /(£") ~ (£" — (/)(0))~". From these it follows 
that the adiabatic growth of a central BH induces a 
density cusp 



p(r) ~ 



S(i?) 

2-7 ' 
4-7 



i?i- 



(14) 



For galaxy models with analytic cores, n = and we 
recover the result A = 3/2, but for models with n > 
we find steeper cusps. We have verified this prediction 
for the five models listed in Table 1. Equation (14) is 
not valid for 7 = 2 (because the potential and distri- 
bution function do not behave as power laws near the 
center), but as 7 approaches 2 the cusp exponent A 
for the 7 models approaches 5/2, which agrees with 
the numerical results shown in Figure 5. 

There is a gap in Table 1 between the models with 
analytic cores, for which A = 3/2, and the 7 models, 
for which A > 2. There is another gap between the 
7 models with < 7 < 2, for which (Dehnen 1993) 



6 



7 



2(2-7)^ 



(15) 



and the model with 7 = 0, for which n = 1. The 
gaps can be filled by other non-analytic models, such 
as the one-parameter family described by the density 



P\(r) 



Po 



(r^ + 



j)4/A' 



(16) 



These models are awkward to work with because the 
potential and distribution function must be found by 
numerical integration. We experimented with some 
models with 1 < A < 2 and found final density cusps 
intermediate between A = 3/2 and A = 2. Equa- 
tion (14) works in some cases but not all; it fails for 
the models with A close to (but less than) 2. 



3.3. Velocity cusps 

The cusps in velocity dispersion in Figures 2-5 all 
rise as {v'^) ~ R~^ at small radii, and do not vary 
much from one model to another, although at low 
BH masses the cusp for the isochrone model is less 
noticeable than for the 7 models. Note that the ve- 
locity dispersion for the 7=1 model without a BH 
(in fact, for any 7 model with 1 < 7 < 2) goes to zero 
at the center (see Binney 1980, Tremaine et al. 1994 
for a discussion of this). 

The anisotropy parameter (3 behaves differently for 
the isochrone model than for the 7 models. The 
isochrone model remains isotropic at the center and 
develops a mild tangential anisotropy away from the 
center, similar to the result for the isothermal sphere 
shown in Figure 1. This is true also for the Plum- 
mer model and, we suspect, for all models with an- 
alytic cores. The 7 models develop larger tangential 
anisotropies that penetrate right to the center. Note 
also that increasing the mass of the central BH has a 
different effect for the two classes of models: for the 
isochrone model it increases the maximum anisotropy 
that develops; for the 7 models it does not change the 
maximum anisotropy, but merely shifts outward the 
radius to which the anisotropy reaches. 

The results for the kurtosis of the LOSVD are more 
difficult to interpret. The figures show the deviation 
of the kurtosis from k = 3, the expected value for a 
Gaussian distribution. In most cases this deviation is 
small in the final cusp, comparable in magnitude with 
what it was in the outer parts of the initial model 
without the BH. Note, however, the following differ- 
ences between the behaviour of the kurtosis near the 
center of the isochrone model and the three 7 models: 
for the isochrone model without a BH the kurtosis is 
constant, whereas for the 7 = 1 model (and other 
7 models with 1 < 7 < 2) the kurtosis diverges; for 
the isochrone model the addition of a BH causes the 
kurtosis to increase, whereas for the three 7 models 
the opposite occurs. 

Perhaps one conclusion to draw from the kurto- 
sis plots is that the adiabatic growth of a central 
BH in a spherical galaxy does not cause the LOSVD 
to become highly non-Gaussian. But note that this 
conclusion applies to the LOSVD measured at one 
exact radius, i.e., to what could be observed if we 
had infinite resolution. The conclusion changes if the 
LOSVD is averaged over an aperture, because the av- 
erage of (lip) will be weighted more towards the center 



(ij = 0) than the average of {vi), and hence the effec- 
tive kurtosis of the averaged LOSVD will differ from 
our comparison of {v"^) and {v£)'^ at the same radius. 
In fact, if the LOSVD is averaged over an aperture 
that includes the BH the effective kurtosis will be in- 
finite, because of the arbitrarily high velocities pos- 
sible close to the BH (see Bahcall and Wolf 1976). 
van der Marel (1994a) shows that this leads to a pos- 
itive /}4 Gauss-Hermite moment, and stresses that it 
is better to quantify the observations by the Gauss- 
Hermite moments than by the classical moments such 
as (v'i) and {v^) . Our work shows at least that, for 
the models we consider, the non-Gaussian nature of 
the LOSVD observed near the center is expected to 
result almost entirely from the spatial averaging, and 
not from a peculiarity of the intrinsic velocity distri- 
bution. 

4. Discussion 

4.1. Theoretical questions 

4.1.1. Anisotropic initial conditions 

The results presented above are for models that 
start with isotropic velocity distributions. We do not 
view this as a severe limitation. Young (1980) sug- 
gested that "in order to have an effect, the anisotropies 
must be significant inside the radius of infiuence r/ ~ 
GMh /cr'i of the black hole as reckoned in the un- 
perturbed cluster." We have verified the correctness 
of this suggestion for Dejonghe's (1987) anisotropic 
Plummer models (we tried models with g = ±1, which 
have anisotropies j3{r) = 0.5(/r^/(l-|-r^)); even for the 
largest BH we considered (Mbh = O.IM) the cusps 
were nearly identical with that for an isotropic Plum- 
mer model. 

We also experimented with some models derived by 
Cuddeford's (1991) technique (with his a set to 1/2) 
to have a constant tangential anisotropy /3(r) = —1/2. 
For the isochrone and 7 = models the result- 
ing cusps differed from those for the corresponding 
isotropic models, but not by much. That is not sur- 
prising, since we know the exact result for the adia- 
batic cusp that forms around a BH if the initial model 
has a density cusp p(r) ~ r~'* made up entirely of cir- 
cular orbits (this is a simple generalization of Young's 
result, derived in Appendix A): 



p{r) ~ 
C = 3 



S(i?) ~ R 



i-c 



3-7 



(17) 



For the isochrone and 7 = models, this "circular" 
cusp slope is C = 9/4, considerably steeper than the 
"isotropic" cusp slope of A = 3/2 for the isochrone 
model, but not much steeper than the A = 2 for the 
7 = model. For the 7 models with < 7 < 2, the 
values of A and C coincide. Since this is for the most 
extreme tangential anisotropy possible, we conclude 
that a moderate tangential anisotropy in the initial 
model will have little or no effect on the final density 
cusp (although it will have an effect on the kinemat- 
ics). 

The case of a galaxy core with a radial anisotropy is 
more difficult to analyse because of a lack of suitable 
models to test. We tried to derive isochrone and 7 = 
models with a constant anisotropy /3 = 1/2 by Cud- 
deford's (1991) technique (with his a set to —1/2), 
but that was not possible: the distribution functions 
turned out to be negative at large binding energies. 
We do not know how large the radial anisotropy can 
be at the center of a galaxy with a fiat core or a mild 
density cusp, but it appears to be small (O. Gerhard, 
private communication), probably too small to signif- 
icantly change the cusp that forms around the central 
BH in our calculations. 

4.1.2. Ongm of BH mass 

Young (1980) suggested that a central BH could 
grow "by accreting gas, from mass loss by giant stars, 
or by some other process," but he added a BH to 
his galaxy models without removing any mass from 
the surrounding stars. We have done some simple 
experiments to check how sensitive the results are to 
this assumption. 

We first tried removing the BH mass uniformly 
from all the stars, by reducing the the stellar dis- 
tribution function by a constant fraction Mbh I M at 
the same time that we added the BH. This made al- 
most no difference to the cusp around the BH, even 
for Mbh /M as large as 0.1. 

We then tried removing the BH mass from the 
stars of lowest angular momentum, since they are the 
stars that approach closest to the center. We adopted 
the following strategy: at the time the BH is added, 
reduce the distribution function f{E, L) by the loss 
fraction If if L < Li, and leave it unchanged if _L > Li, 
i.e., 



(l-lf)f(E,L) 
f(E,L) 



if L < L, 



otherwise. 



f(E,L)- 
with Li chosen so that the total mass removed equals 



Mbh- The rest of the calculation remains the same. 
We assume that, however the mass is lost, it is lost 
slowly so the action variables of the remaining stars 
are adiabatically conserved. The results of several 
calculations of this type for the 7 = model are shown 
in Figure 6. 

The calculation with If = l.O yields a galaxy with 
a hidden BH, with no observable cusp in the surface 
density or projected velocity dispersion. In fact this 
galaxy model has a hole in the middle in two senses: 
a BH and a hole carved out of the intrinsic density 
p(r). Although this is a contrived model — we have 
assumed that all of the mass with L < Li but none 
with L > Lf gets swallowed by the BH, and have ig- 
nored perturbations such as two-body relaxation and 
triaxial components to the potential that might help 
replenish the density hole — it offers the intriguing 
possibility that real galaxies could contain BHs larger 
than suggested by their surface-density and velocity 
cusps. The calculations with If < I show, however, 
that the mass loss must be highly concentrated to- 
wards the center (If ^ 0.5) for it to have a noticeable 
effect. The velocity-dispersion cusp is less sensitive 
to the effects of mass loss than is the surface-density 
cusp. 

4.1.3. Dynamical stability and related questions 

We are confident that our models are dynamically 
stable, but know of no theorem that proves this for 
non-isotropic models with a central BH. The question 
is especially interesting for models with hidden BHs, 
such as that in Figure 6 with If = I. Another ques- 
tion is by how much our models would differ if the 
BHs grow too fast for the assumption of adiabatic 
invariance to be justified. We hope to answer these 
soon with the help of large N-body experiments. 

4.2. Observational implications 

Our results should be compared with observations 
with some caution. We have started our calcula- 
tions from simple galaxy models that, while sharing 
the essential properties of real galaxies, are not ex- 
pected to match them closely at all radii. We have 
assumed spherical symmetry and have ignored rota- 
tion, whereas many of the galaxies believed to contain 
massive BHs have rotating, disc-like nuclei. We have 
presented results that would be obtained with infi- 
nite resolution, and have ignored seeing corrections 
and other such factors that must be included in real- 



istic models. Despite these limitations, we believe our 
results lead to some important conclusions for the in- 
terpretation of density and velocity cusps in galactic 
nuclei. 

4.2.1. Models without BHs 

We start by asking how steep a cusp can be with- 
out a central BH. Dehen (1993) and Tremaine et al. 
(1994) describe 7 models without BHs for all 7 values 
between and 3. Perhaps some of the models with 
large 7 values can be ruled out. 

One restriction is set by the total energy, E = 
— GM/4a(5 — 27), which diverges as 7 approaches 
5/2. Another is set by two-body relaxation. For any 
7 model with 7 > 0, the relaxation time goes to zero 
at the center. This sets a limit to the minimum ra- 
dius for which the model can accurately represent a 
collisionless system (the radius at which t^ equals the 
age of the system). For example, for the models with 
1 < 7 < 3 we find, for r <^ a, 



tr(r) 



0.065t;^ 

M/m 
In A 



1.5 



(3- 

M{r) 

M 



7)(7 - 1)3/2 

3 -7/2 



GM 



1/2 



(19) 



For 7 values close to 3 the exponent (3 — 7/2)/(3 — 7) 
is large, and a substantial fraction of the mass is then 
at radii where tr is small enough to invalidate the 
assumptions on which the model is based (i.e., the 
mass would have undergone core collapse, possibly 
leading to the formation of a BH; see Quinlan and 
Shapiro 1990). While this argument does not give a 
sharp division between acceptable and unacceptable 
values for 7 (because it depends on the values of M , 
m, and a), 7 = 5/2 seems a conservative choice. 

If we therefore disregard models with 7 > 5/2, we 
find that the steepest cusp possible without a BH is 

S(i?) ~ R-''l\ {vD ~ i?-i/2 (^ = 5/2). (20) 

This model has a steep surface-density cusp, much 
steeper than the adiabatic cusp around a BH in an 
isothermal sphere, but only a gradual velocity cusp, 
half as steep as expected around a BH. A surface- 
density cusp at the center of a galaxy thus provides 
only weak evidence for a BH; a Keplerian rise in the 
velocity dispersion towards the center is the real proof 
(if other sources of dark matter can be ruled out). 



4.2.2. Interpretation of density and velocity cusps 

Recent photometry of elliptical galaxies from the 
Hubble Space Telescope (HST) has shown almost 
none to have analytic cores (Crane et al. 1993; Fer- 
rarese et al. 1994; Kormendy et al. 1994); most have 
surface brightnesses that continue to rise at the small- 
est radii resolved. The old picture of a galaxy having 
a core radius r^ within which the surface brightness is 
nearly constant is giving way to a new picture where, 
for many galaxies, the inner parts can be modelled by 
a double power law such as 



i(R) = 2("^-"'y^h ( — 



n 



1+ (- 

n 



(ai-a2)/6 



(21) 

(from Kormendy et al. 1994, with the notation changed 
to avoid confusion with our use of 7), with the tran- 
sition between the two powers ai and a2 occuring 
at a radius rj called the "break" radius or "bend" 
radius or, sometimes (confusingly), the core radius. 
Paradoxically, the larger elliptical galaxies such as 
M87, the ones we expect to harbor massive BHs, have 
only gradual surface-density cusps (ai ~ 0.0-0.3), 
while the smaller ellipticals have the steepest cusps 
(tti ~ 0.5-1.0). The explanation for this dichotomy 
is not clear; it probably requires different formation 
mechanisms, perhaps involving massive BHs. 

The double power law in equation (21) should not 
be confused with the density law (12) for the 7 mod- 
els. These models (such as the Jaffe and Hernquist 
models) were designed to give reasonable fits to ellip- 
tical galaxies if the scale length a is chosen compa- 
rable with the effective radius fg. (Some of the BH 
masses in our calculations are therefore much larger 
than would be expected in real galaxies.) The break 
radius rj in equation (21) is typically a few arcseconds 
for Virgo-cluster ellipticals, at least 10 times smaller 
than re] and the power-law outside the break radius 
is typically a2 — 1.0-2.0, much less steep than the 
~ R~^ fall off for a 7 model at r > a. The 7 models 
are not fiexible enough to match closely a galaxy with 
a double power-law profile at r c:^ ri,. 

Some galaxies have surface-brightness cusps that 
resemble Young's (1980) model for the adiabatic growth 
of a BH in an isothermal sphere. The best known 
of these is M87, first analysed in this way by Young 
et al. (1978). The HST photometry of M87 shows 
a gradual surface-density cusp, I ~ i?"*^^®, which 
Lauer et al. (1992a) say is consistent with Young's 
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model (with Mbh as large as 3 x 10® Mq) because 
the observations do not penetrate to a small enough 
radius to see the expected asymptotic slope of R~^'''^ . 
We checked this with our program and reached the 
same conclusion; observations with a slightly higher 
resolution should see a steepening in the surface 
brightness if Young's model is correct. M32 is an- 
other galaxy studied in detail (Lauer et al. 1992b) 
that can be fit by Young's model, although the fit is 
not as striking as it is for M87. Crane et al. (1993) 
fit Young-type cusps to a number of elliptical galax- 
ies (and give convenient formulas for doing this), al- 
though they find that a single or double power-law 
often fits just as well, sometimes better (especially 
for galaxies with steep inner cusps). 

These fits to Young's (1980) model are suggestive, 
but cannot be accepted as convincing evidence for 
massive BHs. In many cases it is easy to construct 
models without BHs that fit the data just as well. 
We caution against attaching too much significance 
to Young's i?~^" power-law in interpreting surface- 
brightness cusps; steeper cusps result from the adi- 
abatic growth of BHs in galaxies with non-analytic 
cores, and the observations suggest that these are 
common. The adiabatic-growth scenario will always 
cause a steepening in the surface brightness at small 
radii (except in contrived models such as those in 
Fig. 6), although this is not so noticeable for a model 
like Jaffe's that starts with a steep density cusp. 

Convincing evidence for massive BHs can come 
only from high-resolution spectroscopic observations. 
The velocity-dispersion cusp around a BH is insensi- 
tive to the details of the galaxy model within which 
the BH forms (unlike the surface-density cusp which 
varies from model to model), and is difficult or impos- 
sible to mimic without a BH. Ground-based observa- 
tions of M87 show evidence for such a cusp (van der 
Marel 1994b), suggesting that the original Young et 
al. (1978) BH model is correct, a conclusion strength- 
ened by the discovery of a high- velocity gas disk at the 
galaxy center (Harms et al. 1994). More results like 
these from the refurbished HST are eagerly awaited. 
When combined with surface photometry, they can 
help us assess models for the formation of massive 
BHs and dense galactic nuclei. There is still much 
work to do to refine these models so we can extract 
the most information possible from the observations. 
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A. Adiabatic density cusp around a BH 

In the derivation that follows we drop all inessential 
variables and numerical constants (G, tt, M , etc.) and 
consider only the scaling of various quantities with E 
and r. We use subscripts i and / where appropriate 
to distinguish the initial and final states, and adopt a 
sign convention where E and (f) are both positive. We 
approximate the final potential by the Kepler poten- 
tial around the BH. 

Consider first a model that starts with an isotropic 
core. The distribution function f(E) is related to the 
differential energy distribution N(E)dE (the number 
of stars with energies in the range E to E + dE) by 



N{E) = g(E)f(E), 



(Al) 



where the density of states g(E) is, for the isotropic 
initial model, 

■Jo 

Near the BH in the final model the energy varies with 
radius as Ef ~ 1/r, which allows us to relate the 
density pf(r) to Nj{Ej) and Ni{Ei) by 



Pf{r)^r-'Nf{Ef) 



dE 



dr 



'Ni(Ei) 



dEi 



dE 



The relation between Ei and Ef is easy to derive 
for purely circular orbits and purely radial orbits, for 
which the invariance of the action (angular momen- 
tum in the circular case, radial action in the radial 
case) implies that 



p(4-7)/2(2-7) p-1/2 

^i ^f 



(A4) 



E, ~ i?-(2-7)/(4-7) _ ^(2-7)/(4-7), (^5) 

If we use this relation for all orbits we find from equa- 
tion (A3) that 



We thank S. Faber, O. Gerhard, J. Kormendy, 
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Pf(r) 



-r V4-7 



(A6) 
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A model consisting entirely of circular orbits is 
easier to consider. Assume that the density cusp is 
p ~ ^-7 before the addition of the BH and p ~ r'*" 
after. Conservation of mass implies that 

Pirf dr-i = pfrj drj =^ r^ ~^ ~ ''/~'^- i-^'^) 

Conservation of angular momentum implies that 

n Mi {r) = rfMf{r)ct rf Mbh = 



4-7 

r. '~r;. 



Combining these two results we find 

3-7 



C = 3 



4-7' 



(A8) 
(A9) 
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Fig. 1. — Anisotropy induced by the growth of a 
central black hole in an isothermal sphere, as quan- 
tified by (a) the ratio f(E,Lc(E))/f(E,0), and (b) 
the anisotropy /3 in the velocity dispersion. The re- 
sults are from a numerical calculation with a self- 
consistent potential (solid curve), from the same cal- 
culation with an idealized, non-self-consistent poten- 
tial (dashed curve), and from the approximate solu- 
tion of Goodman and Binney (filled squares) . See text 
for details. 



Fig. 2. — Results from the adiabatic growth of central 
black holes of masses 0.001, 0.003, 0.01, 0.03, and 0.1 
in an isochrone model (the dotted lines show the ini- 
tial model without a black hole): (a) surface density; 
(b) line-of-sight velocity dispersion; (c) anisotropy in 
the velocity dispersion; (d) kurtosis (minus three) of 
the line-of-sight velocity distribution. 



Fig. 3. — Results from the adiabatic growth of central 
black holes in a 7 = model (see Fig. 2). 



Fig. 4. — Results from the adiabatic growth of central 
black holes in a 7 = 1 (Hernquist) model (see Fig. 2). 



Fig. 5. — Results from the adiabatic growth of central 
black holes in a 7 = 2 (Jaffe) model (see Fig. 2). 



Fig. 6. — Adiabatic cusp induced by the growth of 
a central black hole of mass 0.01 in a 7 = model 
(cf. Fig. 3) when the black hole grows at the expense 
of low angular-momentum stars. The fractional mass 
loss is If = 1.0 (bottom curves), 0.9, 0.5, and 0.0 (no 
mass loss; top curves). The dotted curves show the 
initial model without a black hole. 
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